﻿/* Copyright 2012 NARI-TECH */
#include <iostream.h>
#include <math.h>
#include <string.h>
#define N 2
double x1, x10, x2, x20;
// double y1,y10,y2,y20;
double y, y10, y2, y20;
double a, b, c, d;
double B[N][N];
double A1[N][N];  // 系数矩阵的逆矩阵

void matrixinversion(double B1[N][N]) {  // 求逆
  int n = 3;        // 统计节点个数（要比雅克比多一维？）
  int i, j, k;      // 循环变量
  double a[N][2*N];        // 求逆矩阵时用的矩阵
  double w = 0, s = 0, l, m1;     // 求逆矩阵时用的变量
  int n1 = n-1;
  for (i = 0;i < n1;i++) {
  for (j = 0;j < n1;j++) {
  a[i][j] = B1[i][j];
    }
  }

  for (i = 0;i < n1;i++) {
    for (j = n1;j < 2*n1;j++) {
       if ((i+n1) == j)
          a[i][j] = 1;
        else
        a[i][j] = 0;
    }
  }

  for (i = 0;i < n1;i++) {
      w = a[i][i];
      for (j = i;j < 2*n1;j++) {
      a[i][j] = a[i][j]/w;
     }
    for (k = i+1;k < n1;k++) {
        s = a[k][i];
        for (j = i;j < 2*n1;j++) {
            a[k][j] = a[k][j]-s*a[i][j];
      }
    }
  }

  for (i = n1-1;i >= 0;i--) {
     l = a[i][i];
     for (j = i;j < 2*n1;j++) {
         a[i][j] = a[i][j]/l;
    }
    for (k = i-1;k >= 0;k--) {
         m1 = a[k][i];
         for (j = i;j < 2*n1;j++) {
             a[k][j] = a[k][j]-m1*a[i][j];
         }
      }
    }

  for (i = 0;i < n1;i++) {
     for (j = n1;j < 2*n1;j++) {
          A1[i][j-n1] = a[i][j];
    }
  }

  cout << "输出系数矩阵B1逆矩阵A1";
  for (i = 0;i < n1;i++) {
  cout << endl;
  for (j = 0; j < n1; j++ ) {
  cout << A1[i][j] << "   ";
}
cout << endl;
}
}

void calculate() {  // 迭代计算
    x10 = 1;
    x20 = 1;
  int i;
  double error = 0.001;
  y10 = 4*x10*x10+x20*x20+2*x10*x20-x20-2;
  y20 = 2*x10*x10+3*x10*x20*x20*x20-3;
  for (i = 0;i < 100;i++) {
  x1 = x10 - (A1[0][0]*y10+A1[0][1]*y20);  // 矩阵相乘
  x2 = x20 - (A1[1][0]*y10+A1[1][1]*y20);
  y = 4*x1*x1+x2*x2+2*x1*x2-x2-2;
  y2 = 2*x1*x1+3*x1*x2*x2*x2-3;
  if ((fabs(y-y10) < error) && (fabs(y2-y20) < error))
  break;
  y10 = y;
  y20 = y2;
  x10 = x1;
  x20 = x2;
  }
  if (i == 100)
  cout << "达到迭代极限" << endl;
  else
  cout << "方程的近似解为" << endl;
  cout << "x1=" << x1 << "   " << "x2=" << x2 << endl;
  cout << "迭代次数为" << i << endl;
  // cout<<y<<endl<<y2<<endl;
}

void main() {
x10 = 1;
x20 = 1;
a = (8*x10+2*x20);   // 雅克比矩阵
b = 2*x20+2*x10-1;
c = 4*x10+3*x20;
d = 3*x10+2*x20;

B[0][0] = a;
B[0][1] = b;
B[1][0] = c;
B[1][1] = d;

matrixinversion(B);  // 求逆

calculate();  // 迭代计算
}